Thermal Conductivity of Anharmonic Lattices: 
Effective Phonons and Quantum Corrections 

Dahai He^j*] Sahin Buyukdagli^j!] and Bambi Hu^'^ 
^Department of Physics, Centre for Nonlinear Studies, 
and the Beijing- Hong Kong- Singapore Joint Centre 
for Nonlinear and Complex Systems (Hong Kong), 
Hong Kong Baptist University, Kowloon Tong, Hong Kong, China 
"^Department of Physics, University of Houston, Houston, Texas 77204-5005, USA 

(Dated: January 21, 2010) 

Abstract 

We compare two effective phonon theories, which have both been apphed recently to study heat 
conduction in anharmonic lattices. In particular, we study the temperature dependence of the ther- 
mal conductivity of Fermi-Pasta-Ulam /3 model via the Debye formula, showing the equivalence 
of both approaches. The temperature for the minimum of the thermal conductivity and the cor- 
responding scaling behavior are analytically calculated, which agree well with the result obtained 
from non-equilibrium simulations. We also give quantum corrections for the thermal conductivity 
from quantum self-consistent phonon theory. The vanishing behavior at low temperature regime 
and the existence of an umklapp peak are qualitatively consistent with experimental studies. 
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The study of heat conduction is very important from both theoretical and experimental 
point of view. A traditional phenomenological approach to understand the thermal proper- 
ties in solids is the Debye formula given by 

K = ^CkVklk, (1) 
k 

where k is the thermal conductivity, Ck,Vk, h are the specific heat, the phonon group velocity 
and the phonon mean free path of mode k, respectively. In spite of its successfulness for 
qualitative explanation of heat conduction in dielectrics, quantitative predictions are hard 
to make from a microscopic viewpoint. 

Recently, an increasing study of heat conduction in low dimensional Hamiltonian models 
may shed light on its microscopic understanding [l| . However, only a few integrable models 
can be solved rigourously Generally one has to rely on numerical simulations for non- 
integrable models. Thus it would be worthy to revisit the traditional kinetic approach by 
incorporating some microscopic consideration for the low dimensional non-integrable lattice 
systems. 

According to the Debye formula, the thermal transport process in a nonlinear lattice is 
intrinsically relative to its dispersion relation and relaxation of normal modes (phonons). 
The existence of nonlinearity makes the definition of phonon delicate, and it is even harder in 
this case to quantify phonon transport from a first-principle way. To surmount the difficulties 
due to nonlinearity, the concept of "effective phonons" has been recently introduced to study 



heat conduction in dynamical models 



The basic idea consists in incorporating the non- 



linearity into normal modes by renormalizing the harmonic frequency spectrum. In the serial 



studies 



-|5|, the authors apply the so-called effective phonon theory (EPT) to study heat 



conduction within the kinetic framework. In Ref. {g], the authors study heat conduction 
through a lattice consisting of two weakly coupled nonlinear segments via the self-consistent 
phonon theory (SCPT). 

It is thus interesting to compare EPT and SCPT since they have both been applied to 
the field of heat conduction. In the present study, we will give a detailed comparison via 



the Debye formula as done 



in la 



-5|. Our result shows the equivalence of SCPT and EPT 
in a large range of temperature. Considering the failure of the classical description at low 
temperature regime, we also compute quantum corrections to the thermal conductivity by 
extending our study to quantum regime, which gives qualitatively consistent results with 



experimental studies. 

It should be emphasized that both EPT and SCPT cannot give the microscopic definition 
of the relaxation time for phonons of mode k. A traditional perturbative way consists 
in studying the single mode relaxation time based on three- and four- phonon processes j3] ■ 
Great efforts have been devoted to obtain the relaxation time of the heat current correla- 



tions 



9|, which in general might be related to the relaxation time of phonons. However, 
the relation was so far unclear. For an anharmonic chain, a simple but physically appealing 
assumption 

T^^ OC WfcC (2) 

has been proposed in Ref. [4], where ook is the phono frequency of mode k. The dimensionless 
nonlinearity e is defined as the ratio between the average of the anharmonic potential energy 
and the total potential energy: 

One can see that, when e vanishes, r approaches to infinity, leading to the expected diver- 
gence of the thermal conductivity. In the following, we will apply the assumption of Eq. ([2]) 
to Eq. ([T]) in order to study the thermal conduction of an anharmonic chain, which gives 
surprising agreement with non-equilibrium simulations. 



I. EFFECTIVE PHONON APPROACHES 



In this section we will compare the SCPT (see, e.g., Q, Inl) and EPT Q Q. To 
make this concrete, we mainly focus on the temperature dependence of the physical quan- 
tities of the Fermi-Pasta-Ulam /3 (FPU-/3) model, which is a classic example to study heat 
conduction. The Hamiltonian of the FPU-/3 model is given by 

^ = E y + f (^^+1 - ^^)' + 1(^^+1 - ^^)'- (4) 

Within the SCP approximation, the Hamiltonian is approximated by a trial Hamiltonian 

where the effective harmonic coupling constant / is given by Eq. flA.29|) in the appendix 
with the Boltzmann constant ks = 1- Note that / is temperature dependent, which stems 
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from the existence of nonlinearity in the system. The dispersion relation of effective phonons 
corresponding to the effective Hamiltonian ([5]) can be written in the form 

UJk = 2^/fsm^. (6) 

It is then straightforward to calculate the dimensionless nonlinearity e and the specific heat 
C from gaussian averages of SCPT, which yields 

f-K 



f + K' 



(7) 



C^U ^ ^ . (8) 

The derivation of the effective phonon spectrum is based on the generalized equipartition 
theorem, 

^sT^{,.^), (9) 

where the bracket (■)^ stands for a thermal average in the canonical ensemble and qu denotes 
the Fourier transform of the coordinate Xj. The next step consists of an approximative 
transformation of the right hand side of Eq. ([9]) into a more compact form (see [s], [l2| for 
details) , 

kBT = ul{ql)^. (10) 

Specifically, the effective phonon spectrum Uk for the FPU-/3 model reads as 

k 

LJk = ^y/a sin-. (11) 

Here a is given by 

= , X (^'+1 ~ ^»)^)^ - 2ATyi/4(x) 

where Yi/4{x) and 1^3/4(0;) are the modified Bessel functions of the second kind and x = 
K'^/{8XT). The parameters for the nonlinearity e and the specific heat C that follow from 
EPT can be expressed in the form 

_ _ K^[Yi/4{x) - Ys/4{x)] + 2XTYi/,{x) 



X 



(13) 



^^3 J<^Yy^ (14) 
4 16ATFi/4(x) 64A2T2^ Yy^ixy' ^ ' 
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FIG. 1: Effective sound speed Vg = \f] as a function of temperature. At high temperature regime, 
Vs oc T^l^. Here K = \,\ = \. 

In Fig. m we plot the effective sound speed Vg = v7 (and Vs = y/a) as a function of the 
temperature. The high temperature behavior gives Vs oc T^/"^, which was already reported in 
Ref. 0, We also compare e and C calculated from SCPT and EFT in Fig. [2] and Fig. |3l 
respectively. One can notice that SCPT and EFT yield practically the same temperature 
dependence over seven orders of magnitude. 

We will now apply these two effective phonon approaches to calculate the thermal con- 
ductivity from the Debye formula (II]). For the sake of simplicity, we will only present the 
derivation from SCPT. The phonon group velocity within SCPT is given by 

Vk = ^ = v^cos I oc v^. (15) 

According to the assumption for the relaxation time (j2]), the mean free path reads as 

h = VkTk oc e"^ (16) 

Substituting Eq. ( TT^ and Eg. ([TB]1 into Eq. ([T]), the temperature dependence of the thermal 
conductivity can be given by |4| 

KiT) oc (17) 
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FIG. 2: Temperature dependence of the dimensionless nonlinearity e. At low temperature regime 
€ ocT, while at high temperature limit e ~ 1. Here K = 1, A = 1. 




FIG. 3: Temperature dependence of specific heat. C ~ 1 at low temperature regime (harmonic 
limit) and C approaches the lower limit 3/4 at high temperature regime. Here K = 1, X = 1. 
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The calculation of k from EPT follows along similar lines. By replacing in Eq. ( !T7|) the 
specific heat, the group non-linearity and the non-linearity by their counterpart deduced 
from EPT, that is Eq. f ll2p . f ll3p and (IT^ . one obtains an alternative expression for the 
thermal conductivity f[T7|) . Note that we only consider the temperature dependence of k 
here as in Ref. [^]. Nevertheless, it should be emphasized that k is both temperature and 
size dependent, k diverges in the thermodynamic limit due to the divergence of summation 
with respect to k in the acoustic regime, which is not the concern of this study. 

The analytical prediction and the simulation result for the temperature dependence of 
K are compared in Fig. [Hand Fig. [5] for A = 1 and A = 0.1, respectively. One can notice 
that the power law behavior observed in Ref. 0, 1^, i.e., k oc at low temperature 
regime and k oc T^/^ at high temperature regime are exactly reproduced. As a comparison, 
non-equilibrium molecular dynamics simulations were performed by applying Langevin heat 
baths at the two ends of the chain [l| . In order to compute the thermal conductivity, a finite 
temperature difference with 10% deviation from the average temperature was consistently 
used. It is clearly seen that the analytical calculations, rescaled by a constant factor that is 
implicit in the Debye formula, give good agreements with the simulation results in a large 
range of temperature. 

Both Fig. m and Fig. [5] show the existence of a minimum of k at a specific temperature 
Tq. By computing the derivative of the thermal conductivity with respect to T, 



(9k 



= 0, (18) 

To 



one easily obtains the turning point in the simple form 

To = go — . (19) 

Here we neglect the temperature dependence of the classical specific heat C, which has a 
very weak dependence on the temperature. The dimensionless constant Qq is then given by 

^0 = ^ + v^, (20) 

If one takes into account the specific heat given by dH]), the behavior of To oc K'^/X survives 
while go should be replaced approximately by 5.52. The turning point Tq denotes the 
transition of k from decreasing behavior to increment behavior with increasing temperature. 
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FIG. 4: Temperature dependence of thermal conductivity. Here K = 1, X = 1, and system size 
N = 1024 for simulation. Analytical results (solid line and dot line) are rescaled (divided by a 
constant), which is consistent with simulation result in whole temperature range. The simulation 
result shows the temperature for the minimum of k at Tq ~ 5, in agreement with Eq. ()19p . 



Eq. ( lT9l) can be understood by scaling the Hamiltonian (jl]) as done in Ref. [14. Let 

K 

Pi = -^Pi, (21a) 
[K 

Xi = d —Xi, (21b) 

one obtains 

H = —H, (22) 
where the dimensionless Hamiltonian H reads as 

^ = E f + ^(^^+1 - ^^)' + ^(^^+1 - ^^)'- (23) 

Eq. ([21]) leads to 

T={p^) = \{PI) = \t. (24) 
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FIG. 5: Temperature dependence of thermal conductivity. Here K = 1, X = 0.1, and = 1024 
for simulation. Like Fig. HI analytical results are rescaled. The simulation result shows Tq ~ 50, 
in agreement with Eq. ()19p . 

Eq. fl2^ means that the temperature dependence of physical quantities of the FPU-/3 model 
can be rescaled according to the temperature transformation in Eq. ( |2^ . which can be 
verified in Fig. O 

Note that the nonlinearity A dependence of k is similar with the temperature dependence 
shown in Fig. HI Fig. [5] and Fig. |6l A critical nonlinearity Aq corresponding to the minimum 
of K can also be similarly obtained as 

Ao = 90 jr- (25) 

Eq. ( |T9|) and Eq. ( l25l) shows that the temperature T and the nonlinearity A are inversely 
equivalent. 

Finally, it is interesting to note that 

e(To) = e(Ao) = (26) 

which is the well-known golden ratio. Due to the particular property of the golden ratio, 
the ratio of nonlinear potential and harmonic potential {Ei)/ {En) is also the golden ratio. 
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FIG. 6: Scaling of thermal conductivity. Here we simply apply coordinate transformation T = T/10 
to Fig. [5] and then merge it with Fig. [H One can see the overlap of the two figures, which indicates 
the scaling of thermal conductivity. 

II. QUANTUM CORRECTIONS 

From the experimental point of view, the thermal conductivity of real materials should 
vanish at zero temperature. This behaviour is inexistent within the classical description of 
heat conduction as shown above, since the latter leads to a divergent thermal conductivity 
at zero temperature. The disagreement between the classical description and experimental 
curves results from the inability of the classical physics to take into account the freezing 
of phonon modes, which is a pure quantum effect. We would like to emphasize that the 
consideration of the low temperature regime is not of pure theoretical interest. Generally, 
the Debye temperature of many solid systems is indeed well above ambient temperature. To 
name but a few, Td = 783 K in stoichiometric CaBg and To = 960 K for vacancy doped 



EuBg [15(1 • In carbon nanotubes, the Debye temperature may even exceed 1000 — 2000 K |16|. 
One thus should take into account quantum effects, namely the partial thermalization of 
phonon modes even far above the room temperature for these materials. 
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In the following, we will treat the Debye equation for the FPU-/3 model semi-classically 
and show that the vanishing behaviour of k manifests itself if one correctly takes into account 
the Bose- Einstein statistics. At this point, we should recall that the derivation of EPT is 
based on the classical equipartition theorem (see Eq. (Q). It is thus impossible to extent 
the effective phonon approach to the quantum regime, which makes SCPT more adequate 
since the latter can be derived from a purely quantum approach. The derivation of the 
quantum self-consistent phonon theory (QSCPT) is presented in the appendix. To avoid 
the divergence of Eq. ([1]) due to the goldstone mode (fc = 0), a small quadratic on-site 
potential /oxf/2 is included in the Hamiltonian (jl]). Here the constant /o is fixed to its 
lower-boundary (/q = 10^^) so that its decrease below this limit doesn't change the value 
of K. A strong point of QSCPT is that the approach unambiguously provides a mode 
decomposition for the specific heat and the phonon velocity, which allows one to consider in 
a rigorous way the discrete mode summation of the Debye equation ([1]). The specific heat 
per mode Ck is obtained from ( 1A.24I) . The phonon velocity Vk follows from the derivative 



of the pseudo-phonon spectrum (IA.21I) and according to Eq. ( JT6l) . the mean free path is 
given by Ik = r^Vk where t^^ oc tUk- For this model, the phonon frequency that should be 
numerically solved with Eq. flA.19|) is 



c^? = /o + 4(i^ + 3A(5a;^»sin2(^), (27) 

where the integer p G [1, iV] and the phonon mode k = 2pTT/N . The non-linearity parameter 
is given by 

fo{x^) /2 + K{5x^) /2 + \{5x^) /A ^ ' 

Fig. [7] gives the specific heat C = 1/A^^Cpasa function of the temperature. At the low 
temperature regime, C increases linearly with increasing temperature since the high energy 
modes are gradually thermalized. This is the well-known partial thermalization effect. 

The agreement with the rigorous result for the harmonic chain (A = 0) at low temperature 
regime shows that the harmonic quantum fluctuations play the dominant role. The quantum 
effect, as expected, is neghgible above the Debye temperature To = 2\/K = 2 and the system 
behaves like the classical one. 

By collecting all the quantities in the Debye equation ([1]), we finally obtain the tempera- 
ture dependence of the thermal conductivity from the quantum approach, which we compare 
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FIG. 7: Specific heat as a function of temperature for (a) A = 1; (b) A = 0.1. The harmonic 
potential is given by Vq = K5x^ /2 + /ox^/2. K = 1 and /o = 10^'' for aU cases. At low 
temperature regime, C ocT. 

in Fig. |8]with its classical counterpart. When the temperature is high enough, QSCPT re- 
produces the classical behavior as shown in Fig. HJand Fig. |5l where the mean free path and 
the specific heat are nearly constant and the temperature dependence of k, mainly follows 
that of the sound speed Vg oc T^/^. The classical and quantum results agree up to the Debye 
temperature. The deviation of the quantum result takes place around Tp, which is followed 
by a peak at temperature T^ax ~ 1, then the thermal conductivity drops to zero. Tmn. 



here indicates the lower bound of the temperature range over which umklapp processes [17 1 
yield the dominant contribution to thermal conductivity. The existence of this peak results 
from a competition between the increasing behavior of the mean free path and the strongly 
decreasing tendency of the specific heat. The latter property is a consequence of the partial 
thermalization of high frequency phonons, which is a pure quantum effect explaining the 
failure of classical approaches below the Debye temperature. At low temperatures, the mean 
free path and the sound speed are nearly constant in the quantum case. The difference for 
the temperature dependence of k is mainly due to the difference of the specific heat. Our 
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FIG. 8: Temperature dependence of thermal conductivity for (a) A = 1; (b) A = 0.1. K = 1 and 
/o = 10^^ for both cases. 



result based on QSCPT, characterized by the existence of the umklapp peak, is quahtatively 



consistent with the experimental studies (see, e.g., Q, Q-ISOl). Finally, one should note 



that the scaling behavior f lMj) is not applicable for quantum case. 



III. SUMMARY 



In summary, SCPT and EFT are compared in the present study. The basic difference 
between SCFT and EFT lies in the different way to calculate the statistical average of 
physical quantities. Specifically, the canonical average is applied in EFT but Gaussian 
average is applied in SCFT. We show that both SCFT and EFT give the same behavior 
for the thermodynamic quantities of a classical system, which indicates their equivalence. 
It seems that gaussian average is preferable for analytical derivation, especially for models 
with an on-site potential, since the thermal averages should be computed in this case from 
numerical transfer matrix calculations. The additional gain is that SCFT can be extended to 
study the quantum system, as shown in section II. Note that EFT is based on the assumption 
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of the equipartition theorem, whose vahdity in the quantum regime is broken. Thus QSCPT 
offers a simple way to explore for a given Hamiltonian model the low temperature physics 
related to the Bose-Einstein statistics. 

We studied the temperature dependence of the thermal conductivity of the FPU-/3 model 
for both the classical and quantum case. The temperature for the minimum of k(T) was 
determined from the classical SCPT, which clearly shows the scaling behavior of thermal 
conductivity for this non-linear model. We also showed that non-equilibrium molecular 
dynamics simulations are in good agreement with the analytical results. At the low tem- 
perature regime, the semiclassical treatment of the Debye equation reproduced the umklapp 
peak, a well-known characteristic observed in experimental studies. 
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Appendix: Quantum Self-Consistent Phonon Theory 

In this appendix, we present a simple derivation of the first order self-consistent phonon 
theory for ID Hamiltonian systems. The presented approach is similar to the derivation of 



the variational path integral method (see 2l|-|23|). In the path integral representation of 



quantum statistical mechanics, the A^-body partition function in the canonical ensemble can 
be expressed as a path integral over periodic trajectories, that is 

Z?xe-^M, (A.l) 



where 

S = j\T[^±' + U[^]) (A.2) 

is the Euclidean action and the first and second terms in the integral correspond to the total 
kinetic and potential energies. It is well known that an exact evaluation of the A^-body 
path integral (lA.ip for general anharmonic potentials is impossible. The idea of SCPT is 
basically approximating the original Euclidean action with a trial action that allows an exact 
evaluation of the trace (lA.ip . Since we deal in this work exclusively with ID Hamiltonian 
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models with closest neighbor interactions of the form 

N 

H = Y.{y^I + ^(^'^) + W{Sxk)} , (A.3) 

k=l 

where Sxk = Xk — Xk-i, an optimal choice of the trial Hamiltonian is that of a coupled 
harmonic oscillator chain, 

i^o = 5^ jyx^ + + AaXfc + s-a^fca^fc-ij • (A.4) 

The trial parameters Ai, A2 and g are to be deduced by minimizing the right hand-side of 
the Feynman- Jensen inequality, 

J'<To+{H-Ho), (A.5) 
where J-q = —ksT In Zq. The trial partition function given by 

Zo = j i^e-^°M/^ (A.6) 
can be easily computed by first performing a Fourier decomposition of the periodic paths, 

n>0 

where Qn = 27Tn/{Ph) stand for Matsubara frequencies and the measure of integration of 
Eq. f lA.6P has the form 



i^=n^^^^=^- (A.8) 

In the last expression as well as in Eq. ( \A.7\i , index k run over oscillators and n over Fourier 
components. We then substitute the expansion (IA.7P into (IA.4P and obtain the trial action 
by performing the integration over time r, 

So= drHo. (A.9) 
Jo 

The trial action Sq that follows is a quadratic function of the Fourier components Xkn, that 
is 

5*0 = S y^fco + ^2Xk + gxkoXk-10 > 

k=i ^ ^ 

N N 

+ /^^E E {(^1 + "^^") + l^^^r) + ^9 KnXlUn + 4>ri„) } (A.IO) 



k=l n>l 
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The next step consists of the trivial integration in Eq. (1A.6|1 over Fourier components, 
which yields 

^°-nsinh(^^,/2)' (^-^^^ 



for Ai + = and 



N 



otherwise. In the last two expressions which give the zero-th order free energy J-q, we have 
defined the pseudo-phonon frequencies in the form 

mul = \i + 2g-Ag sin (A.13) 

and an additional parameter 

V = -Y^- (A.14) 
Ai + 2g 

The calculation of the first order correction to free energy 

{H-Ho) = ^J i^e-^W 1^ iv{xk) + W{Sxk)] - i/o| (A.15) 

proceeds in a similar way. The usual trick consists in expanding the on-site and inter-site 
potential in Fourier basis, 

V{xk) = j ^V{q)e"^^^, W{Sxk) = j ^W{q)e'^'^\ (A.16) 

Then we evaluate the average potentials per particle {V{xk)) and {V{xk)) in flA.lSP by 
integrating over Fourier components and finally invert the fourier transforms of Eq. (1A.16p . 
Consequently, the average value of the potential energies can be expressed in the form of 
smeared-out potentials, 

VM ^ {v{x,)) = J ^7 = {vis^k)) = J y=^^"*^(^)' 

(A.17) 

where we have defined two parameters and 7^ which correspond to lattice displacement 
and two-point correlation function, that is 
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and 

p ^ 



The gaussian smearing f lA.lTp is a key characteristic of variational methods. 
The first order free energy per particle can be finally expressed as 

F, = F„-Y, {^coth } + VM + 11/, (A.20) 

By minimizing this expression with respect to Up, one obtains the optimal pseudo-phonon 
frequency 

On the other hand, the variation of Eq. flA.20p with respects to rj yields 

dV 

^ = 0. (A.22) 
or] 

Note that = for pair on site potentials (y{x) = V{—x)). Using flA.17p . flA.lSp and 
flA.2ip in Eq. flA.20p . we can further simplify the first order free energy in the form 

F^ = Fo + V,-p^^ + W,-,^^. (A.23) 

The specific heat per mode can be then defined by 



Cp = -T^. (A.24) 



One can easily obtain the classical limit of flA.23p by letting — > 0. The final result is 
where and 7^ in Eq. (IA.23P are replaced by their classical counterpart, 
and 

9 /. x9\ 1 v^4sin2(^) 



For FPU-/3 model g]) 



W{5xi) = ^5x1 + j5xt, V = Q. (A.28) 
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Solving Eq. f lA.27p self-consistently with Eq. ( ]A.2ip . it is easy to get the effective harmonic 
constant 

2 _ K + VK'' + 12XKbT 



(A.29) 
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